function logistic

%  comparison of various DE solvers when solving
%  y' = ry(1-y)  with   y(0) = y0     '

% clear all previous variables and plots
clear *
clf

% set parameters for problem
M=10;
r=10;
t0=0; y0=0.1;
tmax=1;

% calculate exact solution
tt=linspace(t0,tmax,100);
a0=(1-y0)/y0;
for it=1:100
	exact(it)=1/(1+a0*exp(-r*tt(it)));
end;

%  calculate numerical solutions
t=linspace(t0,tmax,M+1);
h=t(2)-t(1);
hr=h*r;

% euler
y=y0;
euler=y0;
for i=2:M+1
	yy=y+hr*y*(1-y);
	euler=[euler, yy];
	y=yy;
end;

% b euler
y=y0;
beuler=y0;
for i=2:M+1
	yy=-0.5*(1-hr-sqrt((1-hr)^2+4*hr*y))/hr;
	beuler=[beuler, yy];
	y=yy;
end;

% trap
y=y0;
trap=y0;
for i=2:M+1
	b=2*hr*y*(1+0.5*hr*(1-y));
	yy=(-1+0.5*hr+sqrt((1-0.5*hr)^2+b))/hr;
	trap=[trap, yy];
	y=yy;
end;

% plot results
plot(tt,exact,'k')
hold on
plot(t,euler,'--ro','MarkerSize',7)
plot(t,beuler,'--ms','MarkerSize',7)
plot(t,trap,'--b*','MarkerSize',7)
legend(' Exact',' Euler',' B Euler',' Trap',4);

% define axes and title used in plot
xlabel('t-axis','FontSize',14,'FontWeight','bold')
ylabel('Solution','FontSize',14,'FontWeight','bold')
title(['Logistic Equation:  M = ',num2str(M)],'FontSize',14,'FontWeight','bold')

% have MATLAB use certain plot options (all are optional)
box on
axis([0 1 0 1.1]);
% Set the fontsize to 14 for the plot
set(gca,'FontSize',14); 
% Set legend font to 14/bold                            		
set(findobj(gcf,'tag','legend'),'FontSize',14,'FontWeight','bold');  

hold off